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Abstract. 

We study the response of a two-dimensional hexagonal packing of massless, rigid, 
frictionless spherical grains due to a vertically downward point force on a single grain 
at the top layer. We use a statistical approach, where each mechanically stable 
configuration of contact forces is equally likely. We show that this problem is equivalent 
to a correlated q-model. We find that the response is double-peaked, where the two 
peaks, sharp and single-grain diameter wide, lie on the two downward lattice directions 
emanating from the point of the application of the external force. For systems of finite 
size, the magnitude of these peaks decreases towards the bottom of the packing, while 
progressively a broader, central maximum appears between the peaks. The response 
behaviour displays a remarkable scaling behaviour with system size N: while the 
response in the bulk of the packing scales as on the boundary it is independent of 
iV, so that in the thermodynamic limit only the peaks on the lattice directions persist. 
This qualitative behaviour is extremely robust, as demonstrated by our simulation 
results with different boundary conditions. We have obtained exact expressions of the 
response and higher correlations for any system size in terms of integers corresponding 
to an underlying discrete structure. 
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1. Introduction 

Static properties of granular packings have been a very active field of research in recent 
years. Granular packings are assemblies of macroscopic particles that interact only 
via mechanical repulsion mediated through physical contacts [1,2]. In contrast with 
continuum solids, forces on individual grains in a granular packing can be directly 
accessed both experimentally [3-5] and numerically [6-9]. These forces are found to 
be organized in highly heterogeneous networks that depend strongly on construction 
history of the packing [10]. Statistical studies, motivated to deal with such history 
dependence and heterogeneity of the forces on individual grains, have identified two 
main global characteristics of static granular packings. First, the distribution for the 
magnitudes of inter-grain forces is very broad, with an exponential decay for large 
force magnitudes and a plateau at small force magnitudes [3,6]. Secondly, the average 
response of the packings to a single vertically downwards external force depends strongly 
on the underlying geometry: in ordered packings, the applied force is mainly transmitted 
along the principal lattice directions emanating from the point of application of the force, 
while in disordered packings there is a single vertical propagation direction [4,5,7]. 

Although a large number of different theoretical models have been proposed to study 
these two global characteristics of granular packings, each of these models successfully 
accounts for at most one of them. For example, the so-called g-modcl [11] is a scalar 
lattice model that describes the fluctuations in force transmission within a granular 
packing at the first approximation. While this model does produce realistic distributions 
for the magnitudes of inter-grain forces (see Ref. [12] for a more detailed discussion), it 
predicts a diffusive response to a single external force in conflict with experiments [13]. 
On the other hand, vectorial extensions of the g-model [13], compatible with a more 
general "stress-only" approach [14], were shown to lead to stochastic wave equations. 
For weak disorder, these equations predict a ray-like propagation of stresses in agreement 
with experiments, but for strong disorder, the corresponding behaviour of the stresses 
is less clear-cut. This has further led to the introduction of an ad-hoc "force chain 
splitting" model [15]. The vectorial extensions of the (/-model and the force-chain 
splitting models do predict realistic response behaviour for granular packings, but the 
relation between microscopic stochasticity and resulting distributions for the magnitudes 
of inter-grain forces within these models is not entirely clear. Finally, in contrast 
with these stochastic approaches, classical elastic theory has been used for years in 
the engineering community [16,17]. The predictions of this theory for the response 
behaviour of granular packings match the experimental results rather successfully [18]. 
However, elasticity theory provides only a macroscopic, average level description, and 
thus it provides no information on the force distribution. Moreover, its derivation from 
the grain-level mechanics still seems to be a signiflcant challenge [16, 17]. 

From this perspective, a unifying approach leading to both realistic fluctuations 
in the individual inter-grain force magnitudes and transmissions of forces in a granular 
packing clearly seems to be necessary. Interestingly, in a different context (namely, that 
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of granular compaction) a basis for such an approach has been laid down by Edwards 
years ago [19,20]. By analogy with conventional statistical mechanics, the fundamental 
hypothesis is to consider all "jammed" configurations equally likely. Although no 
microscopic justification for such an ergodic assumption is available so far, experiments 
and numerical studies of quasi-static granular media support such a thermodynamic 
picture [22]. It thus seems very tempting to extend this hypothesis to the study of forces 
in static granular packings by considering sets of forces belonging to all mechanically 
stable configurations equally likely. 

If this idea of equal probability ensemble is to be applied to study the forces in a 
granular packing, it has to be realized that two levels of randomness are generally present 
in the ensemble of forces for stable granular packings [2] . This stems from the fact that 
forces in a granular packing depend critically on the underlying contact network between 
the grains. Contact networks that are deemed isostatic uniquely determine the forces 
admissible on them, and thus, the apphcation of the equal-probability hypothesis to 
packings with isostatic contact networks amounts to considering each contact network 
equally likely. Such a consideration leads to wave equations in disordered geometry 
whach are in conflict with experiments [21]. However, realistic contact networks are 
generically non-isostatic, and the fact that several force configurations can be compatible 
with a given non-isostatic contact network gives rise to the second level of randomness 
[8,9,23]. 

Instead of applying Edward's hypothesis to both levels of randomness 
simultaneously, a natural first step is to apply the uniform probabihty hypothesis first to 
a fixed contact network, and then possibly average over various contact networks. Such 
a study in a fixed contact network has recently been shown to produce distributions for 
the magnitudes of inter-grain forces that compare very well with experiments [24]. In 
Ref. [25] we briefly showed that the application of the uniform probability hypothesis 
in an ordered geometry also leads to a response to an external force qualitatively in 
agreement with experiments. In this paper, we report the same study in full detail. More 
precisely, we determine the behaviour of the response of a two-dimensional hexagonal 
packing of rigid, frictionless and massless spherical grains placed between two vertical 
walls (see FiglH), due to a vertically downward force Wext applied on a single grain at 
the top layer. We define the response of the packing as (Wij) — (VT/j'*) /VText, where 
W,^j and are the vertical forces transmitted by the (i,j)-th grain to the layer 

below it respectively with and without applied external force W^xt- For massless grains 
W^^j = 0. The angular brackets here denote averaging with equal probability over all 
configurations of mechanically stable non-negative contact forces. 

We find that the problem of equally-probable force configurations in the hexagonal 
geometry is equivalent to a correlated g-model. Using this formulation, simulations for 
a variety of boundary conditions show that the response to an external force applied on 
the top of the packing displays two symmetric peaks lying precisely on the two downward 
lattice directions emanating from the point of application of the force. Moreover, 
the response exhibits remarkable scaling with the system size, implying that in the 
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thermodynamic limit the peaks propagate all the way to the bottom of the packing. 
Surprisingly, average values {Wij) as well as higher correlations can be calculated exactly 
for any system size via this formulation. We show that all these quantities can be 
expressed in terms of integers corresponding to an underlying discrete structure. 

The paper is organized as follows. In Sec. El we examine the equal-probability 
hypothesis in a fixed geometry, and its application to the hexagonal geometry in detail. 
In Sec. El we numerically study the response as well as the distributions of the single 
g's and correlations between them. In Sec. jSlwe detail the full theoretical calculation 
of (Wij) and higher correlations for any system size. We finally end this paper with a 
discussion in Sec. (01 




(a) (b) 

Figure 1. Our model: (a) iV x p array of hcxagonally close-packed rigid frictionless 
spherical grains in two-dimensions. At the top, there is only a single vertically 
downward point force applied on one grain, (b) For massless grains, the vertical forces 
are non-zero only inside a triangle formed by the two lattice directions emanating from 
the point of application of Woxt , so that we can restrict our study to that domain. The 
horizontal forces on boundaries of the triangular domain are the same as those on the 
boundaries in (a). 



2. Uniform Measure on the Force Ensemble 
2.1. Force Ensemble in Generic Packings 

To start with, let us examine more closely the inter-grain forces in a granular packing 
[27]. As pointed out earlier, the grains mutually interact only via (mechanical) repulsive 
contact forces. The contacts are assumed pointlike, and in the present study, we consider 
only frictionless grains, so that each contact force is locally normal to the grain surface. 
Thus, the most fundamental entity entering the description of forces in a granular 
packing is the contact network, i.e., the set of all contact points with directions normal 
to the grain surfaces at respective contacts [20]. 
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Consider such a contact network formed by a packing of P grains with Q contacts. 
For simplicity, we restrict ourselves to two dimensions in the following analysis. At each 
contact A; for 1 < A; < Q, the force is represented by its magnitude Fk along the locally 
normal direction to the grain surface, with the convention that a positive magnitude of 
Fk corresponds to a repulsive force. The forces applied on a given grain must satisfy 
three Newton's equations: two for balancing forces in the x and y directions and one for 
balancing torque. For the system in its entirety, the contact forces can be grouped in a 
column vector F consisting of Q non-negative scalars {Fk}, satisfying S — 3P stability 
equations represented by A • F = Fext- Here, Fext is an jS-dimensional column vector 
representing the external forces and torques on the grains, and A is an 5* x Q matrix 
uniquely specified by the contact network. If the grains in question are disks, as is the 
case in most of theoretical and numerical studies, the torque balance is automatically 
satisfied, and the total number of equations S drops to 2P. Note also that there is 
some freedom in the definition of F: a contact force between a grain and a boundary 
can either be considered as an internal force, i.e., as a part of F; or as an external force, 
in which case it becomes a part of Fext- 

In the above description, if A is invertible — which of course implies Q — S 
— the force configuration allowed on the contact network is unique: such a contact 
network is called isostatic. Otherwise cither there is an extended set £ of allowed force 
configurations, in which case the network is called hypcrstatic, or there is none, the 
packing under consideration is unstable under the imposed external forces. 

Theoretical arguments suggest that perfectly rigid grains gcncrically form isostatic 
packings [26]. Physical grains are of course never infinitely rigid, and the corresponding 
contact networks are hyperstatic. These observations are confirmed by numerical 
studies, which moreover find that the convergence to isostaticity with increasing rigidity 
of the grains is rather slow [23] . To determine the forces uniquely, one then in principle 
has to take into account the "force law" , which relates the deformation of a particle to 
the force applied on it, as well as the construction history of the packing. 

Nevertheless, some macroscopic properties of a granular packing, such as the 
ditribution of force magnitudes and the shape of the average response to a point force, are 
independent of the details of the force law [24]. From that point of view, if one considers 
a packing of grains of large but finite rigidity, the deformations of the grains are small 
with respect to their characteristic sizes, and these deformations could be altered without 
essentially modifying the (hyperstatic) contact network [24]. Experimentally, this can 
be achieved by gently tapping the packing without adding or removing contacts [2]. 
Without having to delve deeper in the microscopic force laws, one can then analyze the 
contact force configurations as a subset of the set 8 of allowed force configurations. By 
analogy to classical statistical mechanics, one can study the statistical properties of the 
set £. A natural starting point is to assume that any point in £ is visited with the 
same frequency, similar to a microcanonical ensemble. In other words, one assigns a 
uniform probability measure on under which all allowed contact force configurations 
are equally likely. We should however keep in mind that a priori there is no clear 
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justification for sucli an ergodic liypotliesis. 

Tlie uniform measure on S can be defined more precisely. Tlie set of all solutions 
of A ■ F = Fext is an affine space, whose dimension is the dimension of the kernel of 
A, and £^ is a subset of that space where > \/k = 1 . . . Q. It can be shown 
that £ is usually a compact polygon [8], so that the uniform measure is well defined. 
For a hyperstatic packing, Ker{A) is non-zero, and a parametrization of S can be 
constructed via the three following steps: (1) one first identifies an orthonormal basis 
{F*^')} (I = 1, . . . , dK = Q — S) that spans the space of Ker{A); (2) one then determines 

a unique solution F*^°) of A ■ F*^°) = Foxt by requiring F*^°).F*^'^ = ior I = 1, . . . , dx', and 

Q-s 

(3) one finally obtains all solutions of A ■ F = F^Kt as F = F^^^ + J2 fi F*^'-*, where fi 

1=1 

are real numbers. The restriction of the /^'s to a set S generating non- negative forces 
is a possible parametrization of S. The uniform measure on S is thus equivalent to the 
uniform measure dfi = Y[ dF^ 5(A • F — Fcxt) ©(-^fc) = 11 dfi on S. 

k I 

2.2. Force Ensemble in the Hexagonal Geometry 




Figure 2. Schematically shown forces on the j-th grain in the i-th layer: (a) i = 1, 
(h) i<N and (c) i = N; fI''^' > Vm. 



We will now apply the general method described in Sec. l2.1l to a hexagonal packing 
of monodisperse, rigid and frictionless disks under the locahzed external force We^t (see. 
Fig. P). We will concentrate on massless particles since considering masses does not 
fundamentally change the qualitative behaviour [25]. Our aim is to calculate the mean of 



V3 
' ~ 2 

to be equa' 



(iJ) 
L 

ly likely. 



(see Fig. 12)) by considering all admissible force configurations 



2.2.1. Force balance on individual grains: To start with, we calculate the forces Fg*'"'^ 
and Fq'^^ for each grain in terms of other forces by using Newton's equations. For the 
top layer, i.e., i = 1 [Fig. I2ta)], 
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while for a grain in the bulk [Fig. I21b)] 



5 



(2) 



and finally for the bottom layer [Fig. Et^c)], 

^piNJ) _^ p{N,j) 



p(iy,J) _ pi 



(3) 



In our study, a vertically downward force H^ext is applied on a single grain jo in the 
top layer, i.e., = VFext^jjo- Equations ((UIHl) then show that this force can propagate 

only inside a triangle formed by the two downward lattice directions emanating from the 
grain jo- Outside the triangle, all non-horizontal forces are zero, and for the horizontal 



ones. 



Since our main interest are the 14^(1 j)'s, this implies that we 
can restrict ourselves to the triangular domain, for which the forces exerted on the 
boundaries are the same as the forces on the boundary of the original system [see Fig. 
db)]. 



2.2.2. Parametrization of the force ensemble: As stated in Sec. 12.11 F contains one 
scalar (force magnitude) for each inter-grain contact. This is a consequence of the 
action-reaction principle, which gives the following identifications for 1 < i < and 
Kj <i: 

We consider two different vertical boundaries: (i) "hard walls" described in Fig. 
and (ii) periodic boundaries in the j direction. While in both cases the contact network 
is clearly hyperstatic, the configurations of the internal forces are slightly different. We 
start with the case of "hard walls". 

For reasons which will become clear later on, we will consider the contact forces 
on the right and the bottom boundaries as internal forces, part of F, while 

the forces on the top and right boundaries will be considered as a part of Fgxt- Simple 
counting then shows that Q = ^N^ + ^N, while the number of equations is A^(A^ + 1), 
implying that dx = . 

Notice from Eqs. (0121) that in this description, the lateral forces F3*'"''' and 
F^^'^^ enter the equations for the non-horizontal forces only via their difference Gij = 
j) _ pi^'j) ^Qj. ^ _ 1 ... jv — 1 and j = 1 . . . p. A natural parametrization of S is 
given by these Gjj's: once the Gjj are fixed, all the other non- horizontal forces can be 
uniquely determined by solving Eqs. dUini) layer by layer from top down. It is easily 
seen that the number of these parameters is indeed d^, and that they correspond to an 
orthonormal basis of KerA. 

Thus, the set S is parametrized by the Gjj's, and the Gjj's are restricted within a 
subset S in order to keep all internal forces non-negative. The non-negativity of F^^''^^ 
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»'ox. < Gij. < and - Fl'-'^ < G.J < Fj'-'' . (5) 



Furthermore, the horizontal forces can be expressed as 



pi^.) ^ pi^^) + , (6) 

i'=i 

where F3*'^'*'s are external forces. If the magnitudes of F3*'^'''s are taken to be larger 
then -^Wexti then it is easy to see that F4*'"''''s remain positive for all values of G^j's 

that satisfy inequality (jSj); otherwise, to have F3*'"'-' > 0, the available range for each 

Gij would have to be further restricted within the bounds defined by the inequality 

In this paper, we concentrate on the case ^3*'^-' = -^Wf.^t where the full range (0) is 

allowed. We will consider the case F^'^^ < -^Wext briefly in Sec. 13.11 

Notice that the positivity conditions for F4*'"'^'s defined via Eq. is the only place 
where the precise values of F3*'^'''s enter the analysis. Clearly, considering F^^'^-^'s as a 
part of F would lead to them being parameters of £ (and thus an unbounded S due to 
the lack of upper bounds of Fj^'^'^'s). We therefore choose Fg^'^'^'s as external forces with 
fixed values -^Wey± to keep S bounded and the uniform measure well-defined. 

With the above convention, the set S is an ''^"^''^ -dimensional polyhedron in the 
Gjj-space defined by inequalities The response of the packing to the external force 
VFext is then defined as 

m,)^^jw,,]\dG,,u (7) 

i.e., the averages (. . .) are defined over an ensemble of force configurations, where each 

force configuration is equally likely, with N = I \{dG,, the normalization constant. 

•'^ k,i 

Most of the previous analysis remains valid for periodic boundaries as well. The 
main difference is that the F3*'^'''s are now internal forces (i.e., a part of F), as 
Fg*'^^ = Fi*'^^ due to the action-reaction principle. The phase space S can again be 
parametrized by Gij for i < N,l < j < N, with the additional constraint that 

^G„ = Vz = l...iV. (8) 

i=i 

Once again, only the Gjj 's physically enter the problem, and the precise values of the 
horizontal forces are immaterial. The horizontal forces are well-defined only if one 
chooses fixed reference values for, say, the F3*'^'''s. If these values are large enough, as 
explained in Eq. ® and therebelow, all the horizontal forces are positive irrespective of 
the Gij values within the bounds defined by inequahty Then the set Sp of allowed 
values of the parameters Gij is the intersection of the hyperplane Sh defined by Eq. (jHJ 
and the polyhedron defined by the inequality (jSj). 
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2.2.3. The q- coordinates: There is an alternative formulation of this problem in terms 
of new variables 

where is the fraction of Wi^j that the (i, j)th particle transmits to the layer below it 
towards the left, i.e., Fg*'-'^ = -^QijWij and Fq'^^ = ^(1 - qij)Wij. The mapping Q 
then reduces inequahty (0) to < qij < 1. Clearly, Wij^ = Wg^t is the external force 
applied on the top layer. For i > 1, the Wi/s in the successive layers are related via 

= (1 - qi^ij-i) Wi-ij-i + qi-i,j Wi.ij (10) 

i.e., Wij is a function of qk,i for k < i. 

From the formulation in terms of the g's it may seem from Eq. ()10|) that in the 
hexagonal geometry of Fig. ^ one simply recovers the g-model [11]. There is however an 
important subtlety to take notice of. In the g-model, the g's corresponding to different 
grains are usually uncorrelated. In our case, the uniform measure on S implies, from 
Eq. ©, that 



JV(iV-l) 
2 



l[dQ^,W.AQ), (11) 

i.e., due to the presence of the Jacobian on the right hand side of Eq. (fTTjl. the uniform 
measure on S translates to a non-uniform measure on the -dimensional unit cube 

formed by the accessible values of the g's. In the case of periodic boundary conditions, 
the constraint (jHl) reduces to 

p «-i 2 

E = E = ^^exf (12) 

7 = 1 1 = 1 V<J 



3. Numerical results 



3.1. Shape and scaling of the response 

We now present the results for the {Wij)^s evaluated numerically by implementing 
detailed balance with respect to the probability measure d^ = jj Y{ dqi,j Wij{q) on the 

^) -dimensional unit cube in the g-space. At each step, a single q is modified, and 
the change is accepted with a Metropolis acceptance ratio. From Eq. (fTUjl . it is easy to 
see that the (Wij) values scale linearly with that of W^xt- Thus, from now on, we set 

W^ext = 1. ' 

Our simulation results for (H^jj) in the case of hard walls with ^3*'^-* = -^PVext 
for = 37 are plotted in Fig. Efa), using the built-in cubic interpolation function 
of Mathematica. Outside the triangle shown in Fig. ^b), (PVj.j) = appears in 
deep indigo; the largest (Wij) value within the triangle appears in dark red; and any 
other non-zero {Wij) value is represented by a linear wavelength scale in between. The 
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Figure 3. Simulation results for {Wij): (a) colour plot for (Wij) with hard walls; 
N = 37. The null vanishing forces appear in deep indigo, while the largest are 
represented in dark red. The other values are represented by a linear wavelength scale 
in between. (b)-(c) Behavior of (Wij) in reduced co-ordinates x and z for different A^- 
values and boundary conditions: (b) scaling of {W{x, z)) with system size for |a;| < z 
at two z values; (c) data collapse for {W{x, z)) \x=z for different N values and boundary 
conditions. See text for further details. 



thin green regions that appear on the outer edge of the triangle are artifacts of the 
interpolation. 

We find for any system size N that the (Wi,j) values display two symmetric peaks 
that lie precisely on the two downward lattice directions emanating from the point of 
application of the force Wext- The width of these peaks is a single grain diameter and 
with depth their magnitudes decrease, while a broader central peak appears. At the 
very bottom layer {i = N) only the central maximum for (Wij) remains. 

We further define x = and respectively for even and odd i, and z = 

[see Fig. P in order to put the vertices of the triangle formed by the locations of non-zero 
(Wij) values on (0,0), (—1/2, 1) and (1/2, 1) VA^. The excellent data collapse indicates 
that the (Wij) values for \x\ < z scale with the inverse system size [Fig. Efb); we however 
show only two z values], while the (Wij) values for |x| = z lie on the same curve for all 
system sizes [Fig. Etc)]. The data suggest that in the thermodynamic limit N — > oo, 
the response field {W{x,z)) scales ~ 1/A^ for \x\ < z, but reaches a non-zero limiting 
value on |x| = z Wz < 1. We thus expect ^im {W{x, z))\^^^^^ > {W{x, z))\^^^^^ Wz < 1; 
or equivalently, a double-peaked response field at all depths z < 1 in the thermodynamic 
limit. 

The introduction of the additional constraint p2|) that correspond to periodic 
boundaries modifies neither the qualitative, nor the scaling properties of the response 
function as can be observed in Fig. intb)-(c). 

Although at all places in this paper we consider each side force ^3*'^^ = -^M^xt 
so that the full range of g- values are allowed for each qij, in this paragraph, we take a 
short digression to discuss what happens when one chooses F3*'^'''s to be smaller than 
-^Wext- At the extreme limit when ^3*'^'' = 0, it is easy to see that the only grains 
that correspond to non-zero Wij-values lie exactly on the boundary. From there on, it 
is intuitively clear that once the values of the F3*'^'''s are pushed higher, the finite range 




x/z x/z 



Figure 4. Responses in reduced coordinates for systems with three different values 
of -^i^g'^ at two different depths, with periodic boundary conditions: (a) at z = 0.75, 
(b) at z 1. 



of allowed g-values would begin to transfer some of the vertical forces into the bulk. A 
higher range of allowed values of the g's, caused by higher values of F3*'^'''s, would thus 
effectively result in smaller peaks on the boundary and correspondingly, bigger fractions 
of the total vertical force within the bulk. 

Fig. m where we plot the response at two different z values [z = 0.75 and z = 1) for 
N = 51 and ^^3*'^^ = M4xt, and H^xt/lOO respectively with periodic boundary 

conditions, clearly supports such an intuitive picture. 



3.2. Linearity of response 
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Figure 5. Simulation results for the response N{Wij) due to two vertically downward 
forces of unit magnitude each apphed on the grains (l,ji) and (1,^2) for iV = 35 and 
p = 100: (a) j2 — ji — 35 (b) j2 — ji — 15. Three different values of i are displayed in 
each case. 



Since the values of the {Wij)^s trivially scale linearly with the magnitude of PVext, 
a natural question is whether the response depends linearly on Fext in general, i.e., 
whether the response to a superposition of external forces is simply the superposition 
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Figure 6. p{qi.j) for (a) j ^ i/2 for N ~ 25 and 50, with a quadratic fit in red (b) for 
TV = 25 at i = 15 for four different j vaiues [inset: seif-simiiarity of p{qij) for iV = 25 
and TV = 50 at (x^z) = (0.12,0.4) and(a;, z) = (0.2,0.8)], (c) for TV = 25 on the left 
boundary j = at different heights i. 

of responses. 

It might appear at first sight that the response is indeed hnear — after all, the heart 
of the problem is the system of linear equations A ■ F = Fext- In fact, with the notations 
of Sec. 12.11 one might argue that since F^^-* is a linear function of Fext, while KerA is 
fixed for a given contact geometry, the averages over the affine space F^^-* + KerA should 
depend linearly on Fext- Notice however that the set S corresponds only to F^, > 0, and 
the shape and volume of S in general depend on Fext in a complex way, and thus the 
response needs not be linear. 

To illustrate this point, we apply two vertically downward forces of unit magnitude 
on the (1, ji)th and (1, j2)th grain for N = 35 and p = 100. The |j2 — ji\ > N case is a 
trivial situation, since in this case the two triangular regions of non-zero Wij^s do not 
overlap with each other, and the system is simply the superposition of two stochastically 
independent subsystems. In other words, the response observed for the | j2 ~ ji\ > ^ 
is simply a superposition of the responses of two individual responses [see Fig. Efa)]. 
For 1^2 Ji I ^ N however, the two triangular regions do overlap and the subsystems 
interact; as shown in Fig. Efb), the fact that there exists a central minimum for the 
response at z = 15 is a clear demonstration that the linearity of the responses does not 
hold. 

3.3. q- distributions and correlations 

We have just seen that the qualitative behaviour of our model is very different from that 
of the g-model. It thus seems reasonable that we study the difference between these two 
models at the level of individual grains. 

Our starting hypothesis of equal probability for all mechanically stable force 
configurations implies that the joint probability distribution for the qi/s is given by 
jfY{ijWij{q), where Af is the normalization constant. This distribution does not 
factorize into a product of terms that depend on single gi,j's, and as a result, the qijs 
associated with different grains are correlated with each other. Below we numerically 
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Figure 7. Correlations between Qi^j^ and q.^ j at various values of i for A'' = 50, (a) 
(V, Jc) = (25, 12) and (b) (z„ je) = (47, 10). 



investigate the properties of the induced probabiUty distributions p{qij) for single Qij's, 
as well as the correlations between different Qij^s throughout the system. 

Since individual VFjj(g)'s are independent of qN,kS for k = 1, . . . , N, each of the 
qN,kS is uncorrelated with all other other g's in the system, and is also uniformly 
distributed within the interval [0,1]. However, higher up in the packing, the q- 
distributions start to show a single maximum. For an odd i and jo = i.e., at 
the centre of the layer, p{qijo) exhibits a single maximum precisely at qij^ = 0.5 due to 
symmetry reasons. Furthermore, p{qijo) is independent of the precise value of i < as 
well as the system size, and it is well-fitted by a quadratic polynomial [see Fig. 6(a)]. 
For a given layer i, the larger \j — jo I is, the more the location of the single maximum of 
p{qij) shifts away (symmetrically) from qij = 0.5: for j — jo > 0, this maximum occurs 
at qij > 0.5 and for j — jo < 0, this maximum occurs at qij < 0.5 [Fig. 6(b)]. Finally, 
exactly on the boundary the maximum of p{qi,j) occurs at g = (left boundary) or at 
g = 1 (right boundary) [Fig. 6(c)]. 

Interestingly, the self-similarity that we observed for the response behaviour also 
seems to be vahd for p{qij) in the bulk [see the inset of Fig. 6(b)]. This stands at a 
stark contrast to the physical behaviour of the g-model, as the self-similarity implies 
that in our model, the forces in a granular packing do not propagate from top down, 
instead they depend on the whole extent of the system. 

In Fig. 7, we show results for the correlations between qi^j^ and qi,j for several 
values of i and j for = 50, {ic,jc) = (25, 12) and {ic,jc) = (47, 10). It appears from 
the plots that the correlations between the g's at different locations are very weak. 

4. Exact Results for (Wij): Summary 

Having presented the simulation results above, from this section onwards we concentrate 
on the exact (theoretical) results. It turns out that in this model not only the (Wij) 
values, but also all the higher moments and correlations of the Wij^s can be expressed 
exactly in terms of integers defined by simple recursion relations. The full calculational 
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details for the exact results for (Wij) will be presented in Sec. The details of the 
procedure are slightly involved, and therefore we provide a layout summary of results 
and methods in this section. 

Our main result are exact expressions for all moments of Wij for any system size 
N in terms of integers defined by recursive relations. For instance, the zero-th moment 
or normalization constant J\f for a packing of N layers reads 



{N-1)N 
2 



iN-2) 
■jl,...,JjV-2 



(13) 



where A 



(N~2) 
■jlvj'jV-2 



are integers indexed by {jfc}i<fc<Ar_2, with 1 < ji < 2 and 



^ ^ jk ^ k + jk~i- They are given by the following reccurrence relations 



A 



(2) 
jlj2 



5i. 



a(p+i) 



(14) 



1=1 



A 



(7V-2) 

jl, — ,jN-2 



with 9 the discrete Heaviside (step) function. 

Any higher moment of Wij can be expressed as ;^ I]ji,...,j]v_2 ^j!^...jAr_2' "^^ere 
are integers obtained through slight modifications of the reccurence formulas 
(HD), described in full detail in sections 15.3.21 15.3.31 and 15.41 These results establish 
a non-trivial equivalence between the model studied here and a discrete combinatorial 
problem defined in Sec J5.(il 

The crux of calculating any moment of Wij lies in the transformation of the 
integration measure YidGij [such as in Eq. ((Zj)]. In Eq. (fTT|) . YidGij has been 

rewritten in terms of the Qij^s. With another change of variables from qi j's to VFj^'s, 
the integration measure can be further expressed as 



Y[dG 



JV(iV-l) 
2 



UdqijWijiq) 



JV(iV-l) 



(15) 



While such a variable change simplifies the integrands, the mapped volume S (see 
the last paragraph of Sec. 12.1|1 is however not an A^-dimensional square anymore, but 
again a polygone. Despite this complication, one can build up a recursive structure for 
the integrations over Wk,i 's, which allows to calculate recursively the integrals over S. 

Before we start to calculate any of the integrals, for convenience we first distort 
the triangle in Fig. HJ^b) to that of Fig. IHl^a). We choose jo = 1; and relabel the grains 
in the i-th row as {i,j) with 1 < j < i- Starting at {k,l) = {N,N), we then execute 
the integrals over W^^i in the decreasing order of both k and /; i.e., for a given value of 
k, we integrate over W^^i sequentially for decreasing /, we then decrease k by one, and 
continue the integration process over W^-i,/ sequentially for decreasing /, and so on. 
The integration results over the successive layers are expressed in a "hierarchical nested 
matrix form" in a recursive manner. 
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Figure 8. (a) Schematic view of the sublattice we consider; (b) illustration of the 
"nested" notation used for the matrix elements for p = 2. 



To illustrate the recursive formulation of the "hierarchical nested matrix form", 
let us consider the ealculation of Af. In the expre^ion oi M = f ^<^^^'M. evaluating 

k,l 

the integrations over Wn,n-i through Wn,i, and Wn~i,n-2 through Wn-i,i, we find 
that the result is a polynomial in W^Ar_2,fc that can be written as a matrix product 

„ N-2 N-1 ' ^ N-2 

Jn-2 = / n d^N^i^ill dWN,k = [^^^^] n t^'^'^R^^''- Here, L^^) and R^^^ are 

1=1 k=i 1=1 

respectively 1x2 and 2x1 matrices, and the elements of the 2x2 matrix 

depend only on Wn-2,i- Thereafter, when Jn-2 is integrated over Wn-2,n-3 through 

Wn-1,1, the corresponding result Jn-3 = / II dWN-2,1 Jn-2 yields again a polynomial 

•' 1=1 

r 1 T 

expressible in a similar matrix product form, i.e., Jn-3 = -^'•^^ Yl^^ t^^^'^ R^^\ The 
crucial point to note however is that the matrices L^'^\ t^^^'' and i?*-^-* can be constructed 
by simply unfolding each respective element of L^^\ t^^^'^ and R^^^ as matrices [see Fig. 
IHfb) for the elements of which are now indexed by iiji, 'i2j2]- After the integration 

of Jn-3 over the Vl^Ar_3 /'s, each of the elements of L^'^\ and R^'^^ further unfold 

into matrices and so on. This process continues over further and further integrations 
generating the "hierarchical nested matrix form". The elements of L^'^^^\ j.N-k+1,1 ^^^^ 
flik-i) Q^j-Q related to those of L^''\ t^~'^'^ and R'^^'^ in a recursive manner. 

In Sec. El we develop the full details of this recursive integration scheme and 
calculate {Wij) exactly. In Sec. 15.11 we establish some priliminary formulas, which we 
use over and over again during the course of the exact calculation. In Sees. 15.21 and 15.31 
respectively, we evaluate M and / Hc/py^ ^ How to proceed with the calculations 

k,l 

of the higher moments and correlations of VTij's is discussed in Sec. 15.41 
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5. Exact Results for {Wi,j): Details 
5.1. Preliminaries 

As described in Eqs. f CTTT|) . a full description of a given realization of the contact 
force configurations is given by the variables {qi,j}i=i...N-i,j=i...i with < qij < 1. The 
unnormalized joint probability distribution on these g-variables is Hill^ n}=i ^ij- Here, 
Wi^i = Wext = 1, and for z > 1, the successive Wi/s are given by 

Wij = (1 - qi-i,j-i) W^i-ij-i + Qi-ij for 1 < j <i and 

W^,i = (1 - gi_i,,_.i)iy,_i,,_i . (16) 

We calculate the moments of Wij{q) on this joint g-distribution by changing 
variables from {qij}i-i...N-i,j=i...i to {Wij}i=2...N,j=...i-i- This is achieved by rewriting 
Eq. (Uni) as 



k=l 



+l,fc 



J-1 

fc=i 



(17) 



Notice from Fig. 1(b) that since all the forces on the triangle from the left and the 
right are horizontal, Y^\=iWi^k = IVz. This implies that on the i-th layer there are 
only [i — 1) unconstrained W^jj's. We choose them to be Wij for 1 < j < i — 1. Then 

The advantage associated with the change of variables from {qij}i=i^...,N-i,j=i,...,i 
to {Wi^,} 



jn=2,...,N,j=i, 



.1 is that a=i'n;=i 



nf=Y nr=i w,A ' nf=2 n}=i dw,,^ 



N 



■31 



SO that [from Eq. (fTTjl ] YlijdGi^j = HijdWij, i.e., the Wij's are uniformly distributed. 
The difficulty of this formulation, however, is that the volume the Wij^s span is not a 
cube anymore. Instead, the volume spanned is a polygon given by the inequalities 

(18) 



with 



W, 



i.k) 



i-1 

k=l 

W,_i,,- + ^(iy,_i„ 

k=l 

The ajj's and the bi/s are related by the following simple relations: 



bij 



bi,i 



bij-i - Wij-i 
0. 



Finally, the following integral is the centrepiece of all our calculations: 

rb 



Imn{a,b)= / x"'{b-xTdx = Y.aZ-'' 
•''^ k=0 



ip 



^n+l+fe ^m-k 



(19) 

(20) 
(21) 
(22) 
(23) 

(24) 
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where 



(n + A; + l)!(m- A;)! ' ^ ^ 



5.2. Evaluation of the Normalization constant 

In the Wij-space, the average of a quantity h is given by {h) = J Yiij dWijh /J\f, 
where M = JUijdWij is the normahzation constant. In Sec. 15. 2t we evaluate M by 
performing the integrations over YlijWij layer by layer bottom up; i.e., we start at 
i = N and decrease i until we reach the top layer where i = 1. At each layer, the 
integrations are carried out one by one in the direction of decreasing j, from j = i — 1 
to j = 1. The integral Jn-p = I Y[iLN-p+iT\.j=i,...,i-i dWij is the unnormalized induced 
probability density of the Wij^s for the (A^ — p) top layers. In this notation, M = Ji. 

Our exact calculations are made possible due to the particular forms of the bounds 
ajj's and &ij's, as they introduce a recursive structure. We will evaluate Jn-i and 
Jn-2 explicitly, after which we will prove the general recurrence relation for Jn-p by 
induction. In fact, we will show that Jn-p can be written as a matrix product. The 
matrices entering this product are given by recursive relations, which we will call the 
"fundamental relations" , as they are the building blocks of the calculation of all moments 

„ N-l 

5.2.1. Evaluation of Jn-i CLnd Jn-2- By definition, Jn-i = / II ^Wnj- We 
integrate the Wnj^s one by one in the direction of decreasing j's. Using 



/ dWN,j-i = bN,j-i - aN,j-i = Wn-i,n-i , (26) 



we have 



N-l N-l 



jn-1 = / n d^N,j = n wn-i,j . 



In order to obtain Jn-2, we have to integrate Jn-i w.r.t. Wn-ij for 1 < j < N — 2. 
Since the bounds Uij and bij depend only on Wi^i and VFj-i,/ for / < j, Jn-2 can be 
expressed in a nested form: 



J 



N-2 



N-2 

H dWN-ij Jn-1 

f'JV-1,1 fb N-l, N-2 



/■OjV-1,1 ro N-l, N-2 

/ dWN-l,lWN-l,l ■ ■ ■ dWN-l,N-2WN-l,N-2WN-l,N-l, (27) 

•'ajV-l,! ■^O.N-l.N-2 



and thereafter it can be evaluated iteratively. Having defined ^ = 



N-1,N-1 



1 - T,k=i WN-i,k, we have 



J""-''' = /'"""^ dWN-i,k WN-i,k J"""'''^' (28) 

1jV-l,fc 



Response of a Hexagonal Granular Packing 



18 



for A; > 1, i.e., J, 



J 



N-2A 



. We now show that Vfc > 1, J 



N-2,k 



is of the form 



7l ' (OAr_i,fc_i - WN~l,k-l> + 72 



N-l.k-l 



N-l,k-l 

7i 

N-l,k-l 
72 



and 7 



7 



Af-1.A;-1 



N-l,k 
N-l,k _ 7l 

N-l,k 
72 

^N-l,k ^N-l,k 



where the vectors 7^ ^''^ ^ - 
are related by means of a matrix relation of the form 



For k = N — 1, the proposed form of '^'^ is easily checked by using 
Y.k=o WN-2,k = 1 and Eq. (UHl): 

Af-2 Af-2 

1 - ^ WN-l,k = J2 [^N-2,k - WN-l,k] = bN~l,N-2 " Wn-1,N-2 
k=l k=l 



J 



N-2M-1 



N-l,N-2 r, 1 I Af-l,iV-2 

7i [on^i,n-2 - Wn-i,n^2\ + 12 5 (29) 



with 'ji ^'^ ^ = 1 and 7^ ^'^ ^ = 0. This means that 7^ ^''^ ^ = 

Let us now assume that the form j^-2,A;+i _ j^"^'^ (^b^^i^^ — W^^i^k) + 12 
holds for a given k, 1 < k < N — 2. Then from Eqs. (j24j) and (|28|) . we have 

jN-2,k _ N~l,k J / T \ I N-l,k j / U \ 

J — 7l i'\\\flN-\,k-,ON~\,k) ^ I2 i-\ii\flN-\,k-,ON-\.k) 

0'N-l,k 



N-l,N-2 



1 





N-l,k 1/7 \2 I A^-l.A: 1 /r, \ 

7i ttiil07V-i,fc - aAr-i,fcj +72 aiolOAf-i,fc - CtAf-l,fcj 



+ 



7i 



iV-l,fc 0/7, \3 I N-l,k /T, \5 



Thereafter, using Eqs. ()20j) and ()2ip. ^■'^ can be rewritten as 



J 



N-2,k 



N-l,k-l 



{bN-l,k-l - WN-l,k-l) + 72 



iV-l,fc-l 



with 



N-l,k-l 
7l 

N-l,k~l 
72 



7l v"iv — i.ft— i " iv — i.ft— I /2 5 

s 1 Ti/2 I N-l,k 1 

«llW^7V-2,fc + 72 aioW^Af-2,fc 



N-l,k 1 

7i 



,Af-l,fc Ti/3 



7i *"^"ii'^iv-2,fc + 72 
In other words, we have the matrix relation 
2x2 matrix with elements 

,N-l,k 



'^10^^N-2,k ■ 



(30) 
(31) 

(32) 



^N~i,k^N-i,k^ where t^'^''' is a 



al,2-ii W^7V-2,fc 



(33) 

Thus, by means of the induction procedure via Eqs. ()3UII33|) . we have demonstrated 
that = Uk=i t^~^'^ R^^^ ■ Finally, in the last integral (over Wn-i,i) in Eq. ^ 

the lower limit aAr-i^i = 0, and it is easily seen that 



J, 



N~2 



6]V- 



O-N-1,1 



dWN-l,lWN-l,lJ 



N-2,2 



N-2 



J 



N~2A 



L« n t'^~'''^^'\(34) 



fc=l 



where 



L(^) is the row vector [0 1]. 



5.2.2. Expression of Jjq-p by induction: We now derive the general formula for Jn-p 
by induction. The postulate is that the induced probability density of the N — p top 
layers can be written as 

L^P-W n (35) 



N-p 



k=l 



Response of a Hexagonal Granular Packing 19 
where expressed in "hierarchical nested matrix form", has elements 

,N-ip-l),k _ n ^^rp+ip-i-jp-i 

''iiji,i2j2,..;ip-ijp-i /^nji,«2j2,.--,«p-Up-i '^^ N-p,k K'^^J 

with 

1=2 

Here Q{k) is the discrete Heavyside Step- function, i.e., 9(fc) = 1 if A; > 0, else Q{k) = 0. 
An alternative recursive formulation for the /9's is given by 

Piiji,...,ipjp = Aiii,...,ip_iip_i '^p+jp_i-jp_i,p+jp_i-jp0(v ~ Jp-i) > 1- (38) 

The expressions ()36II38|) will constantly be referred to in all the following calculations, 
and we call them the "fundamental relations" . 

The vectors L^p~^^ and R^p~^\ respectively, are 

p-i 



-^if I 1 = TT ^- 1 1 and 



Mu-}p-i ~ ^jp-iA X! Aiii,.-,«p-2ip-2-RiL--jp-2 • (39) 



/=1 

_ X. - ^ ■ ■ ■ ^^^"^^ 

ii---ip-2 

"Hierarchical nested matrix form" means that elements are referenced through 
nested blocks. For example, the element referenced by ^2^2, • • • , V-iip-i 
ip-ijp-i-th sub-block of the block iiji, 1232, • • • , ip-2jp-2, which itself is the (ip_ijp_i)-th 
sub-block of the block iiji, 22^2, • • • , ip-sjp-s, and so on [see Fig. |Hfb) for an illustration 
in the case p = 2]. The indices vary between the following bounds: 

1 < ?i < 2, 1 < ji < 2 

1<«2<2 + ^1, 1<J2<2+J1 

l<ii<l + ii-i, l<ji<l+ ji-i 



1 < ip~i < P - 1 + ip-2, 1 < jp-i < P - 1 + jp-2 ■ (40) 

In our induction procedure, we assume the above forms (j35H4(jp for a given p, and then 
show that these relations also hold when p is replaced by p + 1 . The two step route for 
this induction procedure that we follow is the same as the one we followed to evaluate 
Jn-2- 

Step (i): We first write 

■^(p-i)]^ (41) 
fc=i 

which we will evaluate in an iterative manner as we did in Sec. 15.2.11 We define the 
vector j^-p-i.^-p = fN-{p-i),N~pj^{p-i) for 1 < A; < — p, and the successive integrals 



J, 



N- 



N-p-1 



n 



N-p,k' 
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are iteratively given by 



jN-p-i,k ^ I dWj,.,,k t^-^P-'^'' jN-p-i,k+i . (42) 



The integral sign in Eq. ()42p is interpreted in the sense that we integrate each component 
of the vector integrand. The final integral (over W^^p^i) yields, just as we saw in Sec. 

Step (ii): We then show that the (ii, «2, • • • , V-i)-th component of J^-p~1'^ is given 

by 

Jt'^ = ' <:t~' (^^-..^-1 - WN-p,k-ir'^-^-'r (43) 

ip=l 

The vectors ^^-P'^-^ and ^^~P'^ with elements and 7^^"^^^ respectively are 

related to each other via the matrix relation = t^-v^^ ryN-p,k ^ 
Starting with j^~p-'^^^-p ^ ^ve have 

jN-p-l,N~P _ \^ n y-r^P+jp_l-fcp_l pP-1 

-'jl,...,jp_i — Pjiki,j2k2,-,jp-ikp-iyv N-p,N-p ^ki,k2,-,kp-i 

ki,k2,...,kp—i 

~ yjiki,j2k2,...,jp-ikp-iyV N-p,N~p'^kp^i,l-n'ki,k2,...,kp-2 

ki,k2,...,kp-i 

= ' 7j:jr;.!!-^ (6^-,,iv-p™i - , (44) 

where we have used Eq. (jH^ between the first and the second lines of Eq. (HH), 
and the fact that Wn-p,n-p = bN-p,N-p-i — Wn-p,n-p-i between the second and 
the third lines. In other words, j^-p-^>'^-p indeed has the postulated form (^^1 with 

= ^jp^ Sfci,...,fcp-i (^hki,...jp-ikp.,RkC!kp_2^ where we choose jp to vary between 
1 and p + jp-i, as shown in Eq. pUjl . 

Assuming the form ()43|) for j^^p-'^'^^ now calculate j^-p-i.'^'-i using Eq. 
The (ii, . . . ,'ip_i)-th element of the corresponding integrand is 

E Au....,v-u.-i<-M7^"W^:!t'(^^-P.^^ - W^.p,k-,r^^-^-^^ , (45) 



so that after using Eq. (|24j) . we obtain 

Jh,...,ip-i ~ E /^niiv,«p-iip-i7ji,..^jp -^p+ip_i-jp_i,p+jp_i-jp('3'Ar-p,fc-l, &Af-p,A;-l)- (46) 

ii,---Jp 

Using Eqs. dSHI) and (|21, we now expand Jp^_jp_-^_jp_-^ p4_jp_j_jp(ajv-p,A;-i) ^Af-p,fc--i) as 

-^p+ip_i-jp_i,p+ip_i-jp(aAf-p,A:-l, i>N-p,k-l) 
p+ip-i-jp-1 

E P+ip-i-jp-i-l P+ip-i-Jp-i-' Ti/P+l+ip-i+'-ip 

"p+ip„i-jp_i,p+jp_i-ip "Af-p,fc-l *'''Ar-p-l,fc-l 

i=0 

p+«p-l 

~ '^l^P Jp~lJ"p+ip_i-jrp_i,p+ip_i-jp "Af-p.fc-l '^'^Af-p-l,fe-l > l^'J 
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where in the last hne of Eq. (j47j) . ip = jp_i + /. Thereafter, we rewrite Eq. (|46p using 
Eq. (jUj) as 



iV-p-l,fc-l 



P+ip~i 



ii,...,i„_i 



X! f^hji,-,ip-ijp-i1ji,.^,'jp X! ®(^P ip-l) ^ 



ir, = l 



^"p+ip_i-jp_i,p+jp_i-ip"Ar-p,fc-l ^^I^Ar-p-l.fc-l 



p+ip-l-ip ^rrP+l+ip-jp 



P+ip-l 

= E7. 



N-p,k-2 p+ip-i-ip 
ii,...,ip ^N-p,k 



P+ip~i 



■ip \^ N—p,k—2 /I -rjT N 

1 — 2^lii,...,ip [0N-p,k-2 - yyN-p,k-2, 



P+jp-i-jp 



(48) 



to recover 



and 



7ii,...,i 



Af-p,fc-2 

,...,ip ~ ''iiji,...,«pjp lji,...,]p 

jl,---,jp 



N-p,k~l N-p,k-l 
7i 



(49) 



,N-p,k-l ^ /I p+ip-i-ip 
iljl,-,ipjp A^njiv,*p-Up-l"p+jp 



1 -ip_ 1, p+ip-l-ip '-'I Jp-lJ ''*'Af-p-l,fc-l 



/=2 



Notice that t^-p-'^-i indeed has the form postulated by Eq. ()36p . 

This procedure outlined in Eqs. ()45II5()|1 allows us to evaluate all the integrals in 
Eq. (PT|) up to = 1. Each iteration in Eq. (PT|) introduces an extra multiplicative 
factor of the matrix t^-p^'^^ whose form is given in Eq. (j50|) . 

For the last integration (over Wn-p,i) in Eq. (jUj), aj^^p^i = and b^-p,! = 1, so 
that simply equal to lfj^;li^_,,p+i^_,; i.e., 



J, 



N-p-l 



«l,«2v«p-l 



ErP-l r ^-P,l 



21,l2,.--i«p 



(p-i) 

il,i2,.--,ip-l"*p.P+«p-l 



*l,«2,---i*p 



AT-p-l 

n t 

k=l 



N-p,k^N-p,N-p 



.N-p-l 



lip) ]J tN-p,k-^ip) 



where 



and 



i.e., Eqs. 



(p) 



k=l 



P 

Lil...,ip_,Si^,p+ip-i = Yl^.^jji+ii^^ 
1=1 



l\,l2,...,lp-\,lp 

(51) 
(52) 



R 



(p) 



N-p,N-p 



A;i,A;2,.--,fcp-i 



d(p-i) 

iifci,i2fc2,.--,ip-ifcp-i-'^fci,...,fcp_i 1 



(53) 



are exactly of the form postulated in Eq. (fHIIjl . 
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5.2.3. The normalization constant Af: The above induction procedure allows us to 
evaluate Jn-p all the way to p = N — 1, and the expression for the normalization 
constant is then given by 



(54) 



For the full expression of Af however, we also need to use Eqs. (|^. and to 
obtain 

^2,1 



nji,«2j2v,«jv-2i]v-2 ~ /^nii>«2i2v,«]v-2ijv-2 since Wi^i — Wext 



(55) 



(2-ji)!(2+ji-j2)! {p+Jp-i-JpV- 



p-i 



{2-ii)\ (2+ii-i2)! "' {p+ip-i-ip)\ {p+l+ip-jp)\ 
and further Eqs. and (jHU)) to derive 



ne(z,-j,„i),(56) 



R 



ip) 

ii,...,ip 



A, 



(p) 



(2-ii)!(2+zi-i2)! . . . ip-l+ip-2-ip-i)l{p-l+tp- 



In Eq. fl57|) . A^^^ are integers given by the following recursive formula: 

"J2,l 

P 



J1J2 



A 



(p+i) 

jl,--;jp+ 



1=1 



Finally, using Eqs, 

X = Ji 



and (jS2I), we have 



{N-1)N 
2 



E a;; 



(7V-2) 



■JJV-2 ■ 



• JlvJiV~2 



(57) 



(58) 



(59) 



5.3. Calculation of the Expectation Values {W^^g^r) 

5.3.1. Modifications of fundamental relations: In this section, we obtain the general 
form of (VF7v-(},r) , the mean value of W^-g^r at q layers from the bottom and r grains from 
the left boundary [see Fig. IHlfor illustration], for < g < — 1 and l<r<A^ — g — 1, 
in four steps. From what we learnt in Sec. 15.21 it is clear that in order to calculate 
(Wj^^g^r), we need to integrate Wjy^g^rJN-q from the {N — q)-th layer to all the way 
to the top. Just like we saw in Sec. 15.21 these integrations are equivalent to matrix 
multiplications layer by layer in an iterative manner, but the fundamental relations 
()36II38|) do get slightly modified. 

Step (i): We start with the calculation of the integral of W]s[-q,rJN-q on the {N — q)- 
th layer 



J, 



N~q-1 



{q,r) 



N-q-l 

n dWN-g,l' 
l'=l 



N-q,r 



N-q 

1=1 



(60) 



wherein once again we carry out the integrations iteratively from V = N — q — 1 down 
to /' = 1. Of them, notice that the integrations for V = N — q — 1 down to /' = r + 1 
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proceed exactly as in Eqs. (|43II5U|) . yielding J^~''~-^'^+^(g, r) = J^-i-^<r+i^ 



JN-g-iiq,r) 



U'=l 
r-l 



n dw, 



N-q,V 



N-q,r 



]-J^iV-(,-l),/ 



u=i 



.1=1 



r . 



(61) 



Step (ii): For I' = r, the recursive integration procedure differs from those in Eqs. 
due to the presence of an extra factor of WN-q,r in the integrand, and thus the 
integration over W^-g^r in Eq- requires a slightly different treatment. We first write 
J'^~'^~^'^{q,r) explicitly: 



fN-q-l,r, 



(iN — q,r 



N-q,r 



X 



^lfj:]^\b^-,.r-Wn-q,rr^^^-^-^^ 
= X! Piljl,...,iq-ijq-l'yji,.!^,'j^ ^q+l+iq-i-jq-i,q+jq-i-jqi(^N-q,r,bN-q,r) ■ (62) 

jl,-,jq 



We then use 



q + l+iq-l-jq-l,q+jq-l-Jq 
q+l+iq-l-jq-l 

= E 

1=0 

(3 + l+jq_l 



l-jq{(^N~q,r: bN~q,r) 

q + l+iq-l-jq-l-l 



^. , - - q+Wq-l-jq-l-l TTTq+l+3q~l+l-jq 

"'N-q,r N-q-l,r 



q+l+iq-^l-iq ^q+l+iq^x-iq TI/'J+l+*9~i' 

q-\-jq-\,1+3q-l-jq 



a,r_.J ' (63) 



(where once again we introduce ig = jq_i + /), to get 

q+l+iq-l 



W+l+ig-i- 



(64) 



with 7^-'?'^ = t(g, 



/ 5 



where 



t 



N-q,r 



q+l + iq-jq 



q + l+iq-jq 



I- n3l,-,tq3q N-q-1, 



(65) 



It is worthwhile to note that all indices except iq in Eqs. 



run between the 



bounds defined in (jlOI), while for iq, we have 1 < iq < q + 1 + iq-i- In other words, the 
fundamental relation (jHBj) gets modified, as t^'i'"^ has one more row than t^~'^'^ . Such a 
modification can be thought of as a "defect" in the recursive integration procedure ()43l - 
ISUj) . More precisely, we call Eq. (jHSl) a "defect of type I at (g, r)" and it is characterized 
by the following relation between (3 and [this relation is obtained via the usage of 
Eq. (I2SD]: 



■t'^q3q 



Alii,.-- 



g + 1 + iq^i - j. 



q-l 



^q3q 



(66) 



q+l + ig-i - 

Step (iii): Next we show that one needs yet another recursive form when V = r — 1 
in Eq. ()60p. but thereafter from V = r — 2 down to /' = 1, the recursive integration 
scheme of Eq. (|6(jp returns to its old form 
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oq,r-l _ n VT-^-rjg-l Jq , . 

f^hjl,i2j2,---iqjq y^lh,nj2,-lqjq „ , o , „' ' ' 



The recursive form for /' = r — 1 is easily obtained by taking the matrix product of 
^N-{q-i),r-i J^~'i~^'^i^q^r) and then integrating it w.r.t. W^-q^r-ii yielding 

Jh..lq-1 (5'' — f^hjl,-,iq-ljq-l ^ q+iq~l~ jq~l,q+'>^+j q-l' jqi^ ^ -q,r~l , ^Af-^.r-l) 

jl,--;jq 

= ' CX"' (^^-.--2 - w^.,,-2r^^-'-'^ , (67) 

iq = l 

where the last line of Eq. (jH7|) has been obtained by using Eq. (j^H). Thereafter, using 
Eq. (I2SI), we get = r)7^~«'^ with 

Ulju...,iqjqiQ^ = Piljl,...,iq-ljq-l ^ ih " J Q^l) '^q+i'q-l- jI-I ,q+l+j q-1- jq ^N-q-l,/-l ' (^S) 

In Eq. all indices except jg satisfy Eq. (HOJ, and for jg, we have I < jg < q+l+jg-i, 
i.e., the fundamental relation is once again modified: r) has one more column 

than t-'v-g.r-i^ pg^jj ^]-^jg modification a "defect of type II at {q,r — 1)", which is 
characterized by the following relation between (3 and (3'^'^~^ [once again, obtained via 
the usage of Eq. ipHjl ]: 

q + 2 + ig-jg 

As for the integrations in Eq. (jUUI) for I' = r — 2 down to /' = 1, notice that the 
final expression (j67|) is the same as the earlier expression (j44j) . although the matrices 
relating the 7's have been modified. Thus, the rest of the integrations in Eq. ()60p for 
the (A^ — g)-th layer yield the recurrence relation ()4H|1 . and the fundamental relations 
for the matrices remain the same as in Eq. ()36|) . The final result then simply becomes 

r-2 N-q~l 

J^_,_i = [l^")] l[t''-'^'^'t''-'^''-\q,r)t''-'^'''{q,r) [] (70) 

l'=l l=r+l 

i.e., in the (A^ — g)-th layer, only two of the matrices get modified w.r.t. Eq. ()36|) . 

Step (iv): In the fourth step, we integrate J^_g_i for the remaining N — q — 2 layers 
to the top of the pile one by one. This integration procedure is the same as what has 
been detailed in steps (i)-(iii), so here we only provide a short description of it. 

When we integrate J^-q-i in the (A^ — Q — l)-th layer, it is easily seen [following 
the calculations in steps (i)-(iii) above] that the fundamental relations for the matrices 
remain the same as in Eq. (j36|l for the integrations over WN-q-i,N-q-2 down to 
WN-q-i,r- Two defects, one of type I and another of type II appear respectively at 
locations (A^ — g — 1, r — 1) and at (A^ — g — 1, r — 2), but the fundamental relations (jHUjl 
are recovered for the integrations over WN-q-i,r-3 down to WN~q~i,i- In other words, 
after the integrations of JN-q-i over all the Wjy^g^i/s for j = (A^ — g — 2), . . . , 1, the 
locations of the defects move one grain each towards the left. In fact, this trend of the 
leftward shift of the defects by one grain each time we integrate over all the Wij^s in 
the successive layer upwards continues to hold until the locations of the defects reach 
(and terminate on) the left boundary [see Fig. IH]. 

In summary, effectively, the difference between the calculations of A/" and that of 
{WN-q,r) lies in the fact that for the latter calculation, the fundamental relations between 
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(q,r) 

N++++++++++ 

Figure 9. Propagation of defects in the recurrence relations for the calculation of 

{WN~q,r)- 



the matrices are modified on the grains located at (A^ — q — k,r — k),0<k<r — 1 and 
(A^ — q — k,r — 1 — k), < k < r — 2. The final result is of the form 

P'\q,r)R^''~'\q,r), (71) 

but the explicit form of (WN^q^r) depends on how the modified fundamental relations 
affect L^^~'^\ r), and R^^~'^\q,r), and hence on the values of r. 

5.3.2. Calculation of {WN-q,r) on the boundary, i.e., r = 1: For r = 1, there is only one 
relevant defect, and it is of type I. It affects only L*^^~^^ and t^'^{q, 1), while R^^~'^\q, 1) 
remains the same as R^^~'^\ Let us first consider the case q > 1, for which we have 

^^•^(?,i)=/5a;!...™-.- (72) 

Herein, the fundamental relation for Aiji,...,ijv-2ijv-2(5'' 1) S^^^ modified only at the g-th 
layer to the form ()66|) to yield 

nq,l ^ n g + 1 + Ig-l - j g-l , . 

A-'nil,...,iiV-2jJV-2 A-'«lil,.--,«JV-2jiV-2 „ I 1 I • _• • V'-'/ 

y ~r J- ~r oq—i oq 

In addition, we should also keep in mind that the dimensions of the matrix 1) 
are not the same as those of t^'^, since the index iq for P'^{q, 1) varies between 1 and 
g + 1 + iq-i as opposed to varying between 1 and q + iq-i for t^'^. This implies that 
the maximum value attained by ii for P'^(g, 1) for Z > g is increased by 1 due to the 
presence of the defect of type I at location (iV — g, r), and consequently 

q-l N-2 

= n ^,„i^-ti)+i n ^jjmi^. m 

l' = l ^ l=q ^ 



{WN-g,r) = if 



liN~2) 
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When Eqs. (17211741) are put together along with Eqs. ((SHI) and ^ and (jSHI) in Eq. (I7T|1 . 

we obtain 

1 



[9(9+1) ,2 ■ 1 


A(.^-2.) 
Jl,---JiV-2 




\N(N-l) ,1 
2 


I 



Jl,---JJV-2 

The two cases g = and g = 1, however, have to be considered separately. For 
g = 1, we find that Eq. (f73j) can be generalized by using = 1, yielding 

while for g = 1, it can be seen that {Wn,i) = (W^7v_i,i)/2. 

5.3.3. In the hulk, i.e., r > 1: For r > 1, each quantity in Eq. (f7T|) is affected due to 
the development of the defects. Of these quantities, t^'^(g, r) is affected by a defect of 
type I terminating on the boundary at location (g + r — 1, 1) and a defect of type II at 
location (g + r — 2, 1). Using Eqs. and (jUIJj) . we then get 

nq,r ^ n g + T - 1 + - jq+r~-2 

Piljl,...,iN-2jN-2 Piljl,...,iN-2jN-2 „ i », _|_ „■ „■ ' 

q -\- r -\- tq-\-r-2 ~ 

We also obtain, just as before, 

q+r-2 N-2 

l' = l ^ l = g+r-l 

However, we still need to express R^'^~'^\q,r). Recall from Eq. (j39|) that i?^^"^'' 
depends, through the recursion formula, on ^^-p.^-p-i for 1 < p < iV — 1. Since for the 
calculation of {W]\f_q^r), t^~P'^~P~i forA^ — r — l<p<A^ — 3 are modified due to the 
defects, R^^~'^\q,r) must also differ from R^^~'^\ 

More precisely, in the recursive formalism described above, R^P\q, r) does not differ 
from R^P\q, r) for p < N — r — 2. The value p = N ~ r ~ 1 corresponds to the layer 
number where a vertical line drawn from (g,r) in Fig. |H1 intersects the right edge of the 
triangle. 

Thus, up to p = — r — 2, the usual recursion formula applies so that 

AN-r-l) 

jp(iV-r-l) _ »lv,»JV-r~l (7Q) 

(2-zi)!(2+zi-Z2)!...(iV-r-2+2jv-r-3-^^-r-2)!(iV-r-2+zjv-r-2)!' ^ ' 
but then is modified with respect to by a defect of type I at (g, r), and we 

get 

_A(^-^) (q r 

-'^-^ (2-ii)!(2+ii-i2)!...(A^-r-l- 

where 



n(Af-r-l)^ q + l+iq-l-iq tl-lN-r--^' ' , . 

n,...,«iv-. ^2-^^)\{2+^^-^,)\...{N-r-l+^,,_r-2-^N--r-l)KN-r-l+^N_,_,)V ^ ^ 



ilvJiV-r-l I 

Once again, it is important to remember that the g-th index ig satisfies 1 < ig < 
g + 1 + while all the others satisfy the inequalities (jiUI) . 
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At the next step of the recursion formula for R^^~'^\ is modified with respect 

to Z?'^''""^ by two defects: one of type I at (g + 1, r — 1) and one of type II at (g, r — 1) 
and we obtain 

I A(^r'-+i) (q r) 

{2-l^)\{2 + l^-l2)\...{N-r + lN-r-l-lN-ry.{N-r + lN^r)\' ^ ' 

where 

a£-:1.(^^0= E \{Qi^l-n-.)^S:l_M^rl (83) 

with 1 < jq < g + 1 + and 1 < ig+i < g + 2 + i^. Notice that Eq. (|H^ is of the 
same form as Eq. except for the bounds of jq and iq+i. 

The calculations for R^^~''\q,r), 2 < k < r proceed along the same lines as for 
R^^~'^^^\q,r) and finally we obtain 

(2-Zi)!(2-2i-^2)!...(Ar-4 + 2^_5-^Ar_4)!(iV-4 + ^^_4)!' ^ ' 

where 1 < iq+r-2 < Q+r—l+iq+rs, and for all the other indices, we have I < k < l+ii-i. 

Putting everything together, as usual, most of the factorials that appear in the 
expressions of P'^{q,r) and R^'^~'^\q,r) cancel out and we are left with 

ilv,jiV-2 

5.4- Higher moments and correlations 

From the methods we described in Sees. 15.1115.31 it is clear that higher moments 
and correlations can in principle be calculated by following the same procedure. We 
will not provide too many details below, instead we will demonstrate that the higher 
moments and correlations can easily be obtained by keeping track of more "defects" in 
the fundamental relations. 

5.4- 1- The case of {W^_g^^) , s > 1 The modifications of the fundamental relations for 
the integration of W^_g^j. are obtained by generalizing the calculations of section 15.3.11 
To start with, the {ii, . . . ,iq-i)-th element of the matrix integration result 

/ ' dWN-q,rW;,_gX-^'^-'^''J'^-'^-'''^\q, r) reads 

E Aiii>---,«g-ii<3-i7ji,...,j, Iq+s+iq-i-jq-i,q+jq-i-jq{0'N-q,r, ^Af-i},r) 

jl,-,jq 

= " <:ribN-q,r-, - W^.q,-,r^^'^-^-^^ (86) 

iq = 

with 7^-'?.-.'^ = t(g,r, s)^-9.'-^Af-g,r+i^ ^^leie 

TN-q,r I „ ^ r,\ — n r\( A \ ^g+'^+iq-l-iq TI/5+l+*q-i^ 



2 +1 



p-g,r ( \ _ o, , . . p)/ - _ ■ q^^^^q-l-^q u/^+H-'g-J, 

hr,s) Tjrq+Wq-jq 
ji,...,iqjq^^N-q-l,r~l 



^^hji,-,iqjq ^N-q-l,r~V ' ) 
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Once again, all the indices vary between the bounds defined in (jlOI), except for ig, 
which satisfies l<iq<q + s + iq-i- The resulting modification of the fundamental 
relations can be called an s-th order defect of type I, and it is characterized by the 
relation 



{q+S + ig-l - jg-l) . . . (g + 1 + ig^l - Jg-l] 



{q + s + ig-i - ig) . . . {q + 1 + ig-i - ig) 

After multiplying the r.h.s. of Eq. ()86|) with the integral with respect 

to WN-q,r-i yields 

^ 7ji,...Jg Iq+iq-i-jq-i,q+s+jq-i-jq{0-N-q,r-l,b]y-q^r-l) 



Jl,---,Jq 



=0 



Furthermore, using 7^ '^'^ ^'^ = t'^ ^[q,r,s)'~f^ we obtain 

-N-q,r-l,Sf ^ _ n p./- _■ q+iq-i-iq u/9+5+l+*'!-i<7 /'qn\ 

''hjl,-,iqjq\'i' ' J — Hti3i,-,tq-i3q-i^\l'q J q~l ) <->i g+i^_^- j ^q+ s+ j q^^~ j q^^ N -q-l,r~l ■> \y^) 

where 1 < jg < s + g + while all the other indices satisfy Eq.ipni). Thus, the s-th 
order defect of type II is characterized by 

5(<?,r-l,s) {q+S+ Jg-l - jg) . . . (g + 1 + - j g) 

Pn..,....,., Pnn,-M3. + , + ^ + _ ^■^) . . . + 2 + - j,) ' ^^'^ 
As in the case for s = 1, all the other integrals on the (A^ — g)-th layer have the usual 
form ()43|1 . so that there are only two defects per layer to take care of. Moreover, these 
defects propagate exactly in the same as shown in Fig. Efa). 

The explicit expression for (W^^r.g in terms of A's can be directly deduced from 
the recurrence relations for the Vs. 



5.4-2. Correlations of the type {W^_g_^^^_^ . . .W^g^^^^J : It turns out that these 
quantities can be calculated using the results for {W^_g^j.). 

Consider the simplest case {W^_g^^j.^W^_g^^^^) with q2 > qi. If the point (q2,r2) 
does not lie on one of the two lines of defects originating from (gi,ri), it is clear that 
the effects on the recurrence relations of the defects originating from (gi, ri) and (g2, 
do not "interact" with each other. 

The defects do "interact" only if (g2, lies on the line defined by (gi + k,ri — k) 
for < A; < ri — 1. In that case, the defects are of Si-th order of type I on (gi + k,ri — k), 
< A; < ri — r2 — 1, and of type II on (gi + /c, ri — A; — 1), < A; < ri — r2 — 2; and 
then of (si + S2)-th order of type I on (gi + k,ri — k),ri — r2 < k < ri — 1. and type II 
(gi + A;, ri — A; — 1), ri — r2 < A; < ri — 2. 

The modifications of the fundamental relations for all the higher correlations 
{^N-qi,ri ■ ■ ■'^N-q„,.r,n) '^^^ obtained by using these observations. Once again if 
the lines of defects originating from the points (g^, r^) for i = 1, . . . , m do not intersect, 
the individual defects do not "interact". If they do intersect, the defects "interact" as 
described above, implying that all correlations can be expressed in terms of A's. 
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Figure 10. Comparison between exact results and simulations: (a) on the boundary 
r = 1 for TV = 11; (b) in the bulk for 9 = 3 and N = 11. 



5.5. Comparison with simulation results 

To evaluate {WijYs exactly using the exact relations developed above for any system 
size, we need to compute the integers A-^-* and correspondingly the K^^^ .-^. These 
integers are defined recurrently by the relations Although the relations ()58p are 

simple sums, the number of terms necessary to evaluate A-^'' increases as roughly as 
(p!)! for large p, so that for practical purposes, it is difficult to go beyond = 11. 
The comparison between the exact evaluation of the (PVjj)'s with the corresponding 
Monte-Carlo simulation results for A^ = 11 is shown in Fig. EB 

5.6. An equivalent combinatorial problem 

We have not been able to find an explicit expression or an asymptotic formula for 
Aj-^ for large p. However, after some relabeling, a combinatorial interpretation can 
be obtained for A^-^^ ^ and related expressions entering the formulas for the moments 
ofW^.,,. 

This interpretation goes as follows: we consider maps h which associate a positive 
integer h(j^^i^ to the (/c, /)-th site of the portion of the square lattice defined hy 1 < k <p 
and 1 < / < A; [i.e., the triangle in Fig. IHl^a)]. The integers hi^k,i) are constrained by the 
following inequahties: 

^(1,1) = 1, 

1 < < 2, V/ = 2...j9 

1 < Vfc) < A; + \/l = 2...p,k = l...l. (92) 

The inequalities are in fact just a re-expression of Eq. (pUj). Furthermore, Eq. 
can then be rewritten as 
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p-1 

A4.i)-^p.p) = '^'^(p>p)'i E A4-M)'-.'^(p-i) n 0[^(p,«+i) - hp~u)]- (93) 

'*(p-l,l)'---'''(p-l,p-l) 

Having iterated the recurrence relation, an equivalent "explicit" expression for A's 
can be obtained as 

Ait,...,,,, = E n nei^w+i) - , (94) 

where the sum runs over all maps h satisfying Eq. (j92|) with fixed values of /i(p,fc) for 
k = 1, . . . ,p. The products on the right hand side are obviously non-zero only for maps 
h satisfying 

h{k,i) < h(^k+i,i+i) Wk = 1 . . .p - 1,1 = 1, . . . ,k. (95) 

We now introduce a new symbol for J2ji,...,jj^^2 ■^^ji,~-^jN-2' which enters the 

expression of the normalization constant JV. It is then easily seen that Z]\}_2 is simply 
the total number of maps h satisfying inequalities (j92|) and for p = — 2. In a 
similar manner, the quantity J2ji,...,jj^_2 Ajf^~j]v_2('?' '^) (jH^ can be re-expressed in 
terms of maps h{q, r) as Z]{ft(g,r)}[Q' + 1 + ^(Af-r,g-i) — ^(Af-r-i,g-i)] where the maps h{q, r) 
satisfy the inequalities and except on the line {l,k) = {N — r + j,q + j) for 
< J < T — 2, where they satisfy 

l<h,k)<k + l + h(^i,k-i). (96) 
Thus, we finally get 

{WN-g,r) = -JjTN^) 7 -^L^^^llA (^g + 1 + h(^N_r^q-l) - h(^N-r-l,q-l))h{g,r)i (97) 

—^-"2 + 1 

where Zf^_2{q,r) is the total number of maps h{q,r) for p = N — 2, and the angular 
brackets on the r.h.s. denote an average over all maps. Similar expressions can also be 
obtained for all higher moments and correlations. 

Interestingly, the original symmetry of the model is not apparent at all in this 
formulation. Although the final results, once calculated [cf. FigHU], display of course 
the symmetry, it seems difficult to show that the underlying discrete problem is indeed 
symmetric. 

6. Discussion and Conclusions 

In this paper, we have studied the response of a hexagonal packing of rigid, frictionless, 
massless, spherical grains to a single external force at the top of it, by supposing 
all mechanically stable force configurations equally likely. We have shown that this 
problem is equivalent to a correlated g-model. Interestingly, while the conventional 
g-model produces a single-peaked, diffusive response, our model leads to two sharp 
peaks on the boundary, i.e., on the two lattice directions emanating from the point of 
application of the external force. For systems of finite size, the magnitude of these peaks 
decreases towards the bottom of the packing, while progressively a broader, central 
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maximum appears between the peaks. The response function displays a remarkable 
scaling behaviour with system size N: while the response in the bulk of the packing 
scales as on the boundary it is independent of N, so that in the thermodynamic limit 
only the peaks on the lattice directions persist. We have obtained exact expressions of 
the (Wij) values, i.e., of the response, and of higher correlations for any system size 
in terms of integers corresponding to an underlying discrete structure, but we have 
not been able to derive an expression for the scaling limit. The resemblance of the 
discrete structure with plane-partition problems [29] might lead to further progress in 
that direction. 

Qualitatively, the response obtained via the uniform probability hypothesis is thus 
in agreement with experiments. It should however be noted that in experiments, the 
width of the peaks increases linearly with depth, while in our model, the peaks are 
single-grain diameter wide at any depth. Such peak widening, in our opinion, is a 
consequence of inter-grain friction. Indeed, a recent study [28] has shown that friction in 
rectangular packings produces a widening of the peaks around the two lattice directions 
emanating from the point of application of the single external force, with the peak 
width proportional to the square root of depth. However, in Ref. [28], the various 
force configurations have not been sampled uniformly as in our model, but rather the 
sampling was carried out in a fashion similar to the usual g-model, with independent 
random values at each grain. We are now attempting to study the effect of friction 
within the uniform ensemble in both rectangular and hexagonal geometries; but our 
preliminary attempts reveal that such a study is technically difficult as the mapping to 
an analogous g-formulation breaks down. 

A natural question that arises in view of our work is whether the uniform probability 
hypothesis leads to one of three categories that the existing continuum-type models for 
the transmissions of stresses have been classified into (namely, elliptic, hyperbolic and 
parabolic according to the nature of the underlying coarse-grained PDE's) [2]. Given 
the critical importance of the underlying network of contacts, it seems difficult to give 
a general answer to this question. For the hexagonal geometry that we studied, it is 
possible to write a coarse-grained equation for vertical stresses u using the g-formulation; 
indeed, since the correlation length between the g's is rather short, the continuous limit 
is the same is in the usual g-model [13] 



where v{x,z) is the noise resulting from the continuous limit of qij = (1 -|- Vij)/2. 
However, in contrast with the usual g-model, the mean {v{x,z)) is position-dependent 
as found in Sec. 13.31 Such a formulation thus does not seem very useful, as the properties 
of the noise and the self-similarity have to be somehow first obtained from the equal- 
probability ensemble hypothesis. Nevertheless, the study of a tensorial formulation of 
this model is an interesting future direction. 

Finally, apart from the inclusion of friction in the present model, the next important 
future direction is the study of the response of disordered granular packings. In that 
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case, averages would have to be taken first for a fixed contact network, and then over 
different contact networks corresponding to different arrangement geometries of grains. 
The shape of the resulting response function would provide another crucial test for 
the applications of the equal probability hypothesis to the study of forces in granular 
packings. 
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